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Abstract 

The instability of a three-dimensional attachment line boundary layer is 
considered in the nonlinear regime* Using weakly nonlinear theory, it is 
found that, apart from a small interval near the (linear) critical Reynolds 
number, finite amplitude solutions bifurcate subcritically from the upper 
branch of the neutral curve* The time-dependent Navier-Stokes equations for 
the attachment line flow have been solved using a Fourier— Chebyshev spectral 
method and the subcritical instability is found at wavenumbers that correspond 
to the upper branch* Both the theory and the numerical calculations show the 
existence of supercritical finite amplitude (equilibrium) states near the 
lower branch which explains why the observed flow exhibits a preference for 
the lower branch modes* The effect of blowing and suction on nonlinear 
stability of the attachment line boundary layer is also investigated. 
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Introduction 


Our concern is with the weakly nonlinear and fully nonlinear stability of 
a three-dimensional attachment line boundary layer obtained by introducing a 
crossflow into the classical Hiemenz stagnation point boundary layer solution. 
The resulting flow has a constant boundary layer thickness and is in fact an 
exact solution of the Navier— Stokes equations. Thus, it is not necessary for 
us to obtain a self-consistent asymptotic solution of the instability problem 
based on a high Reynolds number approximation. In fact, the flow we consider 
is the first order boundary layer solution corresponding to the flow near the 
leading edge of a swept wing. If the flow over the wing is required to be 
laminar, then it is, of course, essential that the attachment line flow be 
stable so that the problem we consider is of direct relevance to laminar flow 
control. 

The present calculation is an extension into the nonlinear regime of the 
work of Hall, Malik and Poll (1984). Hereafter, we refer to that paper as I, 
and we shall shortly discuss the relevant details of that paper. The linear 
theory given in I was motivated by the experimental investigations of 
Pfenninger and Bacon (1969) and Poll (1979, 1980). These authors measured the 
frequencies of naturally occurring disturbances along the attachment lines of 
the flows over different swept cylinders. It was found that small amplitude 
time-periodic disturbances exist above a certain critical Reynolds number and 
correspond to the lower branch of the neutral curve calculated in I. None of 
the experimental points appeared to correspond to upper branch disturbances, 
however, it is, of course, possible that if the flow were forced by, for 
example, a vibrating ribbon, then such modes might be observed. The first aim 
of the present study is to determine whether a weakly nonlinear stability 
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calculation based on the Stuart-Watson expansion procedure can explain why the 
flow exhibits a preference for lower branch modes. We show that, apart from a 
small interval near the critical Reynolds number, finite amplitude solutions 
bifurcate subcritically from the upper branch. This means that these 
solutions are unstable and therefore would not correspond to an observable 
equilibrium state. However, the existence of these solutions suggests that 
the basic state might be nonlinearly unstable to sufficently large finite 
amplitude disturbances. For this reason we decided to investigate numerically 
the full nonlinear stability equations using a Fourier-Chebyshev expansion to 
represent the spatial structure of the disturbance flowfield. 

In fact, Pfenninger and Bacon found that turbulence wires introduced into 
the attachment region could induce large amplitude disturbances in the 
boundary layer at Reynolds numbers significantly below the linear critical 
point. Thus we use a Fourier-Chebyshev spectral method to simulate finite 
amplitude disturbances at Reynolds numbers not necessarily close to the 
neutral curve. In recent years similar calculations for flows such as plane 
Poisseuille flow have become commonplace, and the reader is referred to, for 
example, the papers by Orszag and Kells (1980) and Moin and Kim (1982). This 
type of calculation follows the time evolution of an initial perturbation 
imposed on the basic flow so that unstable time-periodic equilibrium states of 
the type calculated by Herbert (1977) cannot be found by this approach. 
However, the size of such periodic disturbances can be inferred if required by 
gradually increasing the size of the initial perturbation. 

In order to check the results of our calculations, we shall compare the 
numerical results to those predicted by weakly nonlinear stability theory. In 
particular, we calculate numerically the supercritically bifurcating solutions 
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close to the lower branch and see how the size of the equilibrated disturbance 
compares with that predicted by the Stuart-Watson method. Further checks were 
made by comparing our numerical results for small amplitude disturbances with 
the results of I. We shall show that below the linear critical Reynolds 
number it is possible to induce nonlinearly unstable perturbations by 
appropriate choices of the wavenumber and the initial amplitude of the 
disturbance. Qualitatively we will find that our results are consistent with 
the available experimental results. It is possible that the quantitative 
agreement between theory and experiment which we find in the weakly nonlinear 
regime cannot be reproduced in the fully nonlinear regime because the 
disturbances produced experimentally by Pfenninger and Bacon were necessarily 
three-dimensional. The procedure adopted in the rest of the paper is as 
follows: In Section 2 we formulate the stability equations which govern the 
stability of the three-dimensional boundary layer obtained by introducing a 
crossflow into the classical Hiemenz boundary layer solution. In Section 3 we 
discuss the instability in the weakly nonlinear regime whilst in Section 4 we 
discuss the numerical simulation of large amplitude disturbances. Finally, in 
Section 4 we discuss our results and their practical implications. 


2. Formulation of the Problem 

Let us consider the flow of a viscous fluid of kinematic viscosity v 
over the flat plate defined by y — 0. The velocity of the fluid with respect 
to the Cartesian coordinate system (x, y, z) is (u, v, w) and sufficiently 
far away from the plate 


( 2 . 1 ) 
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whilst at the wall we impose the conditions 


u = w = 0, v = V . 

e 


( 2 . 2 ) 


We define A, R, and k by 


4 - (£i) (2.3a) 

e 

W A 

R = — — (2.3b) 



e 


(2.3c) 


so that A is the thickness of the boundary layer at the wall whilst R and 
k are the Reynolds number and a nondimens ional suction parameter 
respectively. 

It is convenient to diverge from the scalings of I and write 


u = W U(X,Y,Z,t), p = pW 1 P(X,Y,Z,t), (2.4) 

~ 6 ^ 6 

where p is the fluid density whilst (X,Y,Z) = A (x,y,z) and the time 
variable t has been scaled on AW^ . The continuity and momentum equations 
then take the form 

V • TT = 0, (2.6a) 

/v 7 

U + (U • V)U = -V P + i V 2 TT, 


(2.6b) 
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and when Y + » we require that U + — , W + 1 . We therefore choose to seek 

K 

a solution of (2.6) which has the particular structure 


U = (XU, V, W), p = - ±_ + P 

R 


(2.7) 


where U, V, W and P now depend only on Y, Z and t. The equations (2.6) 
then simplify to 


u + V y + w z = 0 


U t + u 2 + VU Y + mi z - -i + i ( Uyy + u zz l 

R R 


V t + W Y + m z - - P ¥ + i (v YY ♦ v zz ) 
w t + ™y + TO z * - p z + * I“yy + "zzl 


which are to be solved subject to 


U = W = 0, V=^-, Y = 0 

K 


( 2 . 8 ) 


(2.9) 


U * R ’ 


W ->• 1, Y > 


In the absence of any disturbance, the basic flow takes the form 


U = i u(Y) , V = i v(Y) , W = w(Y) 


where u, v, and w are determined by solving 
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u + v' = 0 


v"' + V ' 2 - v v" - 1 = 0, 


w" - v w' = 0, 


( 2 . 10 ) 


v'(0) = 0, v(0) = k, v'(°°) = -1, w(0) = 0, w(~) = 1. 


When the disturbance imposed on the flow is not small, we must solve (2.8) 
subject to (2.9) numerically, and this will be discussed in Section 4. When 
the disturbance is sufficiently small for us to use weakly nonlinear stability 
theory based on the Stuart-Watson method, it is convenient for us to write 

u-| + tr, 

V = ^ + V, (2.11) 

K 


w = w + w, 


in which case U, V, W and P satisfy 


U + V Y + W z = 0, 


LU - (2u U + V u'| = R{U 2 + VU + WU }, 

X Lt 


LV - Vv' = RP y + R{VV Y + WV Z }, 


( 2 . 12 ) 


LW + RVw' = RP Z + R{VW Y + WW 2 }. 
Here the operator L is defined by 




(2.13) 


The discussion of I was restricted to the linear regime where the nonlinear 
terms in (2.12) can be neglected, in the following section we determine how 
these disturbances develop in a neighborhood of the neutral curve whilst in 
Section 4 larger disturbances will be calculated by integrating (2.8) 
numerically. 


3. Weakly Nonlinear Stability Theory 

In I the solution of the linearized version of (2.12) was discussed, this 
was done by taking U, V, W and P to be proportional to 

E - exp ia[Z - ct] . Thus the disturbance has wavelength 2ir/a and 

propagates along the attachment line with speed c. We found in I that in the 
case of zero suction instability is possible for R > 583.1 and that with 
suction the flow is significantly stabilized. We follow the usual approach of 
weakly nonlinear stability theory and determine how the disturbance develops 
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in a neighborhood of a point on the neutral curve in the a-R plane. Suppose 
then that (a 0> R Q ) is a point on a neutral curve for some values of k and 
that the corresponding value of c is cq. We expand 

R = R q + eR^ + •••, (3.1) 

where 0 < e « 1 and define a slow time variable t by 

t = et. (3.2) 

The X velocity component then expands as 

U = {e 1/2 UqE + eU 2 E 2 + e 3/2 U 3 E 3 + e 3/2 U 4 E} + COMPLEX CONJUGATE + eU M + 0(e 2 ) 

and V, W and P are expanded in a similar manner. It is then a routine 
procedure to substitute the above expansions into (2.12) and equate like 
powers of e ^ 2 . At order e ^ 2 we find that 

<VV ' A(T)< W 

where A(t) is an amplitude function to be determined at higher order whilst 
(U 0 ,V 0 ) satisfies the sixth-order differential system 
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(Mi + i« 0 K q c 0 }U 0 - 2u0 0 - »'V 0 - VU' - 1« 0 ^ » D 0 - 0 

( M 1 + lo o R 0 c o) M l V 0 ~ la 0 R 0 “ M 1 v 0 + ia 0 E 0 ”" v 0 ~ 7m i v o ■ »' M i v o 

+ 2u'U q + 2uU' + U" V 0 + U'V' = 0 

u 0 = V 0 = V 0 = °» Y = °»“- (3.3) 

_ d 2 2 2 

where j - a Q j . Thus (3.3) is just the eigenvalue problem discussed 

in I, and in Figures 1 and 2 we have shown the neutral values of a a c 

0’ 0 0 

for several different values of k. The eigenrelation was obtained by using 
a fourth-order accurate finite difference scheme to solve (3.3) after first 

writing V = [ V g > v q > v q^ > v o " » U o ,U ()] T so that Z satisfies an equation of the 
form 

dV 

d? =A l (3.4) 

where A is a 6x6 matrix whose elements are given explicitly in I. Later we 
shall need the solution of the system adjoint to (3.3) and if 

3. ” , ^2’^3 , ^4 ,< ^5 ,< ^6^ t * ie adjoint vector, the appropriate system is 

d_q 

dy = “A q, q 3 = q 4 = q 6 = °» Y = 0,». (3.5) 

Here the precise manner in which these functions decay to zero can be found by 

looking at the asymptotic solution of the adjoint differential equation for 
Y » 1. We note here only that if we insist that this decay is exponential, 

then (3.5) has only a discrete spectrum of eigenvalues which, of course, is 

identical to that associated with (3.3). 


r 
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At order e we find that 


(u 2 > v 2 ) - a 2 (b 2 ,v 2 ), (u m> v m ,w m> p m ) - |a | 2 (u m ,v m ,w m ,p m ). 


where (u 9 ,V„) satisfy 


(M 2 + 2ia 0 R 0 c 0 }B 2 - 2uU 2 - u'V 2 - vB' - 2i<. Q R fj wB 2 - R^U' - B Q V'}, 


(m 2 + 2ic 0 R 0 c 0 )M 2 V 2 - 2io 0 R 0 wM 2 V 2 + la Q R Q »"V 2 - vM 2 V' 


■ - R 0 [- 4a o “o v o + 2 ( 2u o "8 + v o V + a 'o v 'o - V 0 V 


- “o “o' + 2u o v o')l’ 


(3.6) 


whilst the mean flow correction is determined by 


“m + “m ■ 0 


(3.7a) 


V - - 2 "“m - v m + (V 0 U 0>' + (V 0 V'l 


(3.7b) 


v » - Vv m - V M - “o - “ol« v o v o>' + “o v o + “o v o) 


(3.7c) 


w m' - vW l 


M - R 0 V M W ' = 1 ST < v 0 <+ v 0 - V 0 u 6 - V o' V 0 >> (3 * 7d) 


where * denotes 'complex conjugate.' The system of equation (3.6) is to be 


solved subject to 
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u 2 = v 2 = v 2 = °» Y = 0, 

with U 2 , V 2 tending to zero exponentially when Y -*■ <». Now turning to the 
mean field correction equations, we note that U M can be eliminated from 
(3.7b) using the equation of continuity to give a third-order equation for 
Vm* For large values of Y this equation has the three independent solutions 

(i) V. . ~ constant, 

M 

(ii) V m ~Y 3 , 

(iii) V M ~ e _Y /2 

so that in order to satisfy the no-slip condition V M = U M = 0 at the wall, 
we relax the condition on V M at 'infinity' to 

V' + 0, Y i- «. (3 . 8) 

Having determined V M , we can then integrate (3.7d) to find W M and we note 

that the structure of the equation for W M for Y » 1 enables us to find 
such that 


W M (0) = 0, 0 (exponentially), Y ®. 


(3.9) 


Finally, and Pj^ can then be determined from (3.7a,c) respectively. 

3/2 

At order e we obtain differential systems for (u^> Vj, Wj> P 3 )» 

/A/ /v /v rss v 

1^4 > ^ 4 » ^ 4 » P 4 J » *-he usual way. We obtain an amplitude equation for A(t) 
as a solvability condition on the system for (u 3 , V 3 , W 3> P 3 ). The equation 
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takes the form 


dr a 0 R 1 A + a l A ^ 


where the constants aQ and aj are defined by 


- f [%(<» - c o ) ( v o' - “o V + "" v oK + - c o ) 0 ok] dY 


R 0 l 0 [ U 0 «6 + K - “0 V oKJ dY 


‘1 i"[ D 0’6 + (V - »0 V 0^4] dY 


‘ - f R o [ r 4 !’ “o( V 0 V i + V M V i + ia 0 W M V 0 + 3V 2 V o' + I V 0 V 2 


+ 2V 2 D 0 + 1 U 2 V o) - ld ( V oV - U M U 8 - W i “0 - Vo ^ 


+ { (2vJ'u' + V* U'' + V* V'" + 2V* V" + 2VJ, U' + 2V„ U" 


+ 2 V' V" + 2V V'""’ - 2V„ U„ 
M 0 M 0 2 0 


’k*’ k' 

- V 2 °0 - 2V 2 V 0 


- v' v 
2 0 


+ U 2 u*' + U' uj + U 2 V* + u* V")j+ q 6 ji«o w M°o + v M u o + v o u t 


+ 2 U 0°M + i U 0 U 2 + V 0 U 2 + 1 U 0 V 2 + D 0 v 2 + 2 U 2 V o’i I"’ 


(3.11b) 
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The constants a Q and a± can be determined only after integrating 
numerically the differential systems for the eigenfunction, adjoint 
eigenfunction, first harmonic function and the mean field correction. This 
was done using the fourth-order accurate finite difference scheme described in 
I, in all the calculations the eigenfunctions were normalized such that the 
maximum value of |Wq| = 1. The integrals appearing in the definition of aQ 
and a 1 were then evaluated using the trapezium rule, the results of a 
calculation for k = -.1, 0, .4, .8 are shown in Figure 3. We have only 
given results for the real parts of aQ and a^ since this information is 
sufficient to calculate the amplitude of an equilibrium disturbance. Before 
discussing these results further we note from (3.10) that 

7 ■ 37 W 2 - R 1 a 0r W 2 + a lr W 4 

so that equilibrium solutions are possible if 


and this solution is stable if the bifurcation is supercritical (aj r < 0) 
and unstable if the bifurcation is subcritical (a lr >0). In the latter case 
a finite amplitude motion having | A | ^ aQ r a^£ causes | A | to increase 

without limit. 

Now let us turn to the results illustrated in Figure 3, the most 
important results correspond to k = 0, and we see that the bifurcation is 
always supercritical on the lower branch. But a^ r has a zero near 
r q ~ 595, and for the remainder of the upper branch the bifurcation is 


7 
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subcritical. In Figure 4 we have shown the neutral curve for tc = 0 together 
with the experimental points due to Poll (1979, 1980) and Pfenninger and Bacon 
(1969). We have marked by an arrow the position along the upper branch beyond 
which the finite amplitude solution bifurcates subcritically and is 
unstable. As expected it seems that the experiments have picked up the 
disturbances corresponding to the parts of the neutral curve where the 
bifurcation is supercritical. In Figure 5 we have shown the eigenfunction, 
first harmonic and mean flow correction corresponding to the critical point on 
the neutral curve for k = 0. hater we shall describe a numerical 
investigation of finite amplitude disturbances, and in order to test our 
calculations we shall try to reproduce quantitatively the finite amplitude 
solution which bifurcates supercritically from the lower branch at R = 800. 
We shall also investigate the possibility of finite amplitude motions at 
Reynolds numbers significantly less than the critical value, these 
disturbances are to be expected since the bifurcation is subcritical over most 
of the upper branch. 

It remains for us to discuss the results foe the cases when k * 0. We 
see that increasing the blowing at the wall reduces the Reynolds number regime 
over which subcritical disturbances are possible. In fact when k = .8 the 
bifurcation is always supercritical so that the flow is not susceptible to 
'threshold amplitude' effects. 

We see in Figure 3 that when k = -.1 the point on the neutral curve 
where there is a crossover from subcritical to supercritical bifurcation moves 
down from the upper branch to a point on the lower branch. Thus, the 
bifurcation is now subcritical at the critical Reynolds number, this is 
consistent with the results of Hocking (1974) who investigated the nonlinear 
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stability of the asymptotic suction boundary layer which corresponds to the 
limit k Furthermore, we note that at some values of k in this range 
(-.1,0) the zero of a lr must occur at the critical Reynolds number; thus, 
if we wanted to extend our analysis to include slow spanwise variations the 
appropriate evolution equation would not in this case be that found by 
Stewartson and Stuart (1971) and would have to include fifth-order terms in 
the disturbance amplitude. 

The evolution of the supercritically bifurcating solution with increasing 
Reynolds number is beyond the scope of the present calculation. However, if 
the disturbances develop in a manner typical of convective or centrifugal 
instabilities, it is possible that the flow remains laminar over a significant 
range of values of the Reynolds number. If three-dimensional instabilities of 
the supercritically bifurcating solution exist then the subsequent development 
of the flow would be more complex. However, if the origin of transition on 
the attachment line of a swept wing is due to subcritical disturbances, then 
it is not clear whether suction should be effective in keeping the attachment 
line stable. This follows from the fact that suction, although increasing the 
critical Reynolds number, makes the flow more susceptible to subcritical 
disturbances. In contrast, blowing ultimately causes the disappearance of 
subcritical disturbances but lowers the critical Reynolds numbers at which 
infinitesimal disturbances are unstable. 


4. Direct Numerical Simulation 

The attachment line boundary layer is strictly parallel; i.e., the basic 
flow is independent of the coordinate along the attachment line. Therefore we 


-16- 


can employ periodic boundary conditions in that direction for the solution of 
(2.8). This is to be contrasted with the Blasius boundary layer where 
periodic boundary conditions do not simulate the actual physical problem in a 
rational way as the growth of the boundary layer cannot be accounted for. 

For the present boundary layer, a Fourier-Chebyshev spectral method will 
be used to simulate two-dimensional finite amplitude states. We use the 
spectral-collocation method of Malik, Zang and Hussaini (1984) (hereafter 
referred to as MZH) for the solution of (2.8) subject to the boundary 
conditions (2.9). A stretching transformation can be applied in the 
(unbounded) vertical direction. Let 


Y - 


+ 

b - n ’ 


(4.1) 


where Y is the physical vertical coordinate, n the computational coordinate 
and a and b are constants. Let Y max be the upper boundary in the 
physical plane and set 


b = 


1 + 


2a 


Y 

max 


(4.2) 


Then for any choice of the scaling parameter a, the computational 
coordinate n falls within the standard Chebyshev interval [-1,1]. The 
collocation points in the computational plane are 


Z j " J V K > 


j = 0,1, •••,K-1 


n 


m 


m=0, 1 , • • • ,N, 


(4.3) 


/■miT'v 


(4.4) 
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where L z = 2ir/ct, and K and N are the number of intervals in the Z and 
Y directions, respectively. The dependent variables have Fourier-Chebyshev 
series of the form 

K/2 N 2ulk/L 7 

u(Z,n,t) = £ E u kn (t)e T n (n) > (4.5) 
k= -K/2 n=0 kn n 

where T n is the Chebyshev polynomial of degree n. In the spectral 
collocation method, spatial derivatives of u are obtained by differentiating 
the series expansion with the expansion coefficients u kn (t) determined by 
discrete Fourier and Chebyshev transforms of the grid point values of u. The 
details of the procedure are given in Gottlieb and Orszag (1977). Derivatives 
in the vertical direction are evaluated by multiplying the Chebyshev 
collocation in n by the Jacobian of the transformation, i.e., 

“Y = "3Y u n* < 4 * 6 ) 

In the temporal discretization, the pressure gradient term and the 
incompressibility constraint are best handled implicitly. So, too, are the 
vertical diffusion terms because of the fine mesh-spacing near the wall. We 
use Crank— Nicholson time discretization on the implicit terms and second— order 
Adams-Bashf orth on the remainder. After a discrete Fourier transform in Z, 
the following set of ordinary differential equations result (we list them in 
the order they are stacked for numerical computations): 
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u n+1 - vf 1 - ikW n+1 


BV ^ 1 + v n+1 + qJ +1 = V n + 


H?" 1 ) 


Qy + ev" Y 


BW ^* 1 + W n+1 + ikQ n+1 


W n + ^ (3H? 


H^" 1 ) " ikQ n + PW YY 


- bu ^ 1 + u n+1 


J n + ^| (3H" 


sr 1 ) + <. 


(4.7) 


A /V A 

In the above, k = 2irk/L z , B = At/ 2R, Q = AtP/2, i = /=T, and “ denotes 
Fourier transformed variables in wavenumber space. The wavenumber is denoted 
by k and the dependence of W, V, Q and U (the order of the dependent 
variables here represents that of the solution vector adopted for the 
numerical solution of (4.7)) upon k has been suppressed. The superscript 
n represents the time level* 9 H 2 and H 3 , which contain the terms 

treated explicitly, are given by 


Hj - -W Y - WV Z + 5 v zz 


H 2 - -vw ? - w z + T w zz 


H 3 - - VU Y - WU Y + i U zz - 0 2 + 


Appropriate boundary conditions are yet to be prescribed for (4.7) which will 
be discussed later. 

For each wavenumber k, the system of equations (4.7) can be written as 
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L S = F, (4.9) 

where S = [w n+1 , V n+1 , Q n+1 , U n+1 ] and F is the known right-hand side. 
The matrix L constructed by using Chebyshev polynomials is a full MxM matrix 
where M = 4N. A direct soltuion of (4.9) by Guass elimination methods would 
require O(M^) storage and O(M^) arithmetic operations. In the present 
study, we use a spectral iteration scheme, based on a minimum residual (MR) 
method with finite difference preconditioning, that requires only 0(M) 
storage and 0(M log M) operations per iteration. An effective 
preconditioning is provided by using a staggered mesh in the normal direction 
whereby the velocities are defined at the cell faces n m , and the pressure at 
the cell centres 


V-l /2 = cos(ir(m- i /2 )/N), m= 1,2, ••*,«. (4.10) 

The momentum equations are enforced at the cell faces, whereas the continuity 
equations are enforced at the centres. More details of the iterative spectral 
method employing staggered mesh are given in MZH. 

We now return to the question of boundary conditions imposed for the 
numerical solution of (4.7). Because of the staggered mesh in the vertical 
direction, no artificial pressure boundary conditions are required. The 

a 

velocity boundary conditions for k * 0 are 

U = v = W = 0, 

and 

^ A A 

u = V = W = 0, 


0 , 


max’ 


(4.11) 
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or 


U Y * - Y «ax U > V Y ‘ - l k l V > W Y ’ * l k l H > 


Y = 


Y 

max 


In some test runs, both the zeroth-order and first-order boundary conditions 
at Y = Y„,_„ gave almost identical results when Y > 15. With the 
first-order boundary conditions, the iteration scheme (MR) converged faster; 
and therefore, these conditions were imposed at ^ max = ^ in t,ie 

calculations to be reported in this section. This fast convergence with 
first-order boundary conditions was also noted in MZH. 

For k = 0, the boundary conditions are 


V 



W = U = 0, 


Y = 0 


W = 1, 



Y = 


Y 

max 


(4.12) 


The structure of (4.7) for k = 0, with the above boundary conditions, is 
quite simple. In this case W and U satisfy two' tridiagonal equations, and 
after first solving this system the continuity equation is then solved as a 


bidiagonal equation for V. Once V is known, the pressure Q also 

A 

satisfies a bidiagonal equation. This is solved by setting Q(Y i j ) =0 and 
then solving for each successive value of the pressure. This particular 

A 

choice of Q(Y y ) is arbitrary and corresponds to specifying the mean 

pressure. 

Initial conditions required for the solution of (4.7) are provided by 
imposing a disturbance of finite amplitude upon the basic state. The 
disturbance eigenfunctions are calculated using linear theory as discussed in 
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I. The initial conditions thus are 


U(Z,Y,0) = ^- + e Re{U 0 (Y)e iaZ } 

V(Z,Y,0) = + e Re{V 0 (Y)e iaZ } (4.13) 

W(Z,Y,0) = w(Y) + e Re{W () (Y)e ictZ } > 


where the disturbance eigenfunctions U Q , V Q , W Q have been normalized such 
that the maximum value is 1. The parameter e has been introduced to control 
the magnitude of initial disturbance. 

Let us define the flow energy at any time t as 


E(t) 



(u(Z,Y,t) -£) 2 + (v(Z,Y,t) -^) 2 + (w(Z,Y,t) - w) 2 JdY 

(4.14) 


and rate of change of the disturbance amplitude as 


a 


2E dt 


(4.15) 


where a > 0 and a < 0 signify growing and decaying disturbances 
respectively. 

In the numerical calculations which are reported below, we have used 33 
Chebyshev polynomials in the normal direction whilst the number of Fourier 
modes along the attachment line varies from case to case. Excellent agreement 
between the numerical results and linear theory was achieved in MZH for plane 
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Poiseuille flow and Blasius boundary layer with only 33 Chebyshev polynomials. 

We compared the solution of (4.7) with linear theory results just to check the 

accuracy of the numerical scheme and found satisfactory agreement. As an 

example, calculations were performed at a = .25648 for three different 

Reynolds numbers using e = .0001 and K = 4. The results are presented in 

L Z 

Table 1. These results were obtained using a time step At = CFL — with a 
CFL number of 0.1. Calculations were terminated at t = 306 which 

corresponds to about 4.5 linear wave periods. For all three Reynolds numbers, 

the difference between the calculated a (averaged) and linear theory result 
is approximately .000016 which is indicative of the degree of accuracy that 
can be expected with the employed spatial and temporal discretizations. In 
order to estimate numerical dispersion in the calculation scheme, we performed 
a computation at R = 570 with a = .32 and e = .00001. The calculated 

wall pressure for this wave is plotted in Figure 6. The nondimens ional 

frequency calculated from the signal is .1235: the corresponding linear 

theory result is .1249. Having established that reasonably accurate results 
may be expected from the numerical computations when N = 32 and CFL = 0.1, 
we now present some results that pertain to finite amplitude motions. 

According to linear stability theory, all infinitesimal disturbances 
decay for k = 0 if R < 583.1. The critical wavenumber in this case is 
a = .288. The weakly nonlinear theory presented in Section 3 showed that 

bifurcation is always supercritical near the lower branch of the neutral curve 
and is subcritical on the upper branch so the flow will be unstable in a 
finite-amplitude sense for wavenumbers corresponding to most of the upper 
branch of the neutral curve. We first show that the numerical computations 
support the result that subcritical bifurcations cannot take place at 
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wavenumbers that correspond to the lower branch of the neutral curve. We do 
this by performing a calculation at R = 570 with a = .28, e = .12 and 
K = 8. The results are presented in Figure 7(a) - 7(d). In Figure 7(a) the 

disturbance energy is plotted and is found to decay as a function of time. 

Figure 7(b) contains a plot of the rate of change of disturbance amplitude 
(o) which shows that after an initial period of positive growth, a settles 
down at a negative value of about -.00059 (the linear value is -.00016). 
Amplitudes of the fundamental mode and first harmonic are plotted in Figure 
7(c) which also decay. Finally, a trace of the calculated wall pressure 

signal is plotted in Figure 7(d) . Similar calculations were performed at R = 
570 with a = .25648 and e = .05, .12, .2. The results are consistent with 
the weakly nonlinear theory result that unstable (subcritical) finite- 

amplitude disturbances cannot exist in a swept attachment line boundary layer 
at wavenumbers that correspond to the lower branch of the neutral curve. Our 
full nonlinear computations do support the prediction of the weakly nonlinear 
theory that subcritical instability can occur at wavenumbers that correspond 
to the upper branch of the neutral curve. Our computations at R = 570 with 
e = .12 and a = .32, .33, .34 and .37 all show the existence of unstable 
finite-amplitude motions. The band of unstable wavenumbers at R = 570 lies 
in the range .28 < a < .4 with a = .34 as the most unstable wavenumber. 
The results for this wavenumber are presented in Figure 8(a) - 8(d). Figure 
8(a) shows that flow energy increases with time. The disturbance growth rate 
is plotted in Figure 8(b). The magnitude of the growth rate at the time when 
computations were terminated was about .00036: the corresponding linear 

theory result is -.00099. The amplitude of the fundamental mode and first 
harmonic are plotted in Figure 8(c) while the wall pressure distribution is 
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given in Figure 8(d). These results were obtained using K = 8; however, some 
of the computations at other wavenumbers were done using K = 16, and very 
little effect on the growth rate was found. These computations clearly show 
the existence of subcritical instability in the attachment line boundary 
layer. In Figure 9 the effect of varying Reynolds number is studied for 
e = .12, o = .34 and K = 8 where the disturbance energy E(t) is plotted 
for R = 570, 550, 540, and 530. It appears that for a = .34 and e = .12, 
subcritical instability is present only for Reynolds numbers R > 540. At 
R = 530, E(t) decays when e = .12. Calculations at other wavenumbers also 
showed a similar trend. Higher initial disturbance amplitude may, however, 
trigger subcritical instability at lower Reynolds numbers. Some calculations 
performed at R = 538 showed that the growth rate increases with increasing 
e. 

A set of calculations was carried out at R = 500 with e = .15 and .2; 
the disturbances at all the wavenumbers decayed. In Figure lOa-d we present 
the results of these calculations for a = .35, K = 16. The growth rate at 
the end of the computation is about -.00051 whilst the corresponding linear 
value is -.00197. A summary of all the computations at subcritical Reynolds 
numbers is presented in Table II. Based on these computations it appears that 
subcritical instability could exist in swept attachment line boundary layer 
only at Reynolds numbers in excess of about 535 with reasonable amplitudes of 
initial disturbances. The experimental points of Pfenninger and Bacon (1969) 
below this Reynolds number may have been influenced by three-dimensional 
disturbances. 

We now come to the question of the supercritical bifurcation near the 
lower branch of the neutral curve as predicted by the weakly nonlinear 
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theory. For the bifurcation to be supercritical, a neutrally stable solution 
should exist at supercritical Reynolds numbers. We investigated this result 
by performing a computation at R = 900 and a = .201178. For these 

conditions the linear growth rate is .00065. We chose an initial amplitude 

e = .04 and a streamwise resolution with K - 8. As shown in Figure 11(a) 

the energy solution settles down at a constant value. This is confirmed by 

plotting the growth rate a (Figure 11(b)) which attains a value of 
approximately zero at large t for e = .04. The amplitude of the 
fundamental mode is plotted in Figure 11(c). The final value is .0398 in 
excellent agreement with the weakly nonlinear prediction of .04. In order to 
see the effect of initial disturbance, two solutions were obtained with 
e = .02 and e = .06 which are also plotted in Figure 11. It is quite clear 
that they too tend to the same finite— amplitude solution. 

Another calculation was performed at the same wavenumber but at a higher 
Reynolds number of 1000. According to the weakly nonlinear theory, the 
equilibrium amplitude should be about .056 for this Reynolds number. First, a 
computation was performed with e = .05. The amplitude of the fundamental 
mode is plotted in Figure 12. At the end of the computation, the amplitude 
is slightly above .06 and still increasing. The calculation was repeated with 
e = .065 and the result is also plotted in Figure 12. We see that this 
solution shows an equilibrium state at an amplitude of about .067. The 
discrepancy between this and the weakly nonlinear prediction is not totally 
unexpected since at this Reynolds number the contribution from the higher 
order nonlinear terms will be significant. 
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5. Conclusion 

The results which we have given in this paper extend the linear theory of 
I into both the weakly and fully nonlinear regions using asymptotic and 
numerical means. There seems little doubt that the absence of any upper 
branch modes in the experiments of Pfenninger, Bacon, and Poll is due to the 
subcritical nature of the bifurcation along most of the upper branch. 
Furthermore, this subcritical bifurcation is the cause of the nonlinear ly 
unstable disturbances which are investigated numerically in Section 4. These 
disturbances exist below the linear critical Reynolds number and the region in 
the wavenumber-Reynolds number plane where they are unstable increases with 
the size of the initial amplitude. The existence of these modes is 
consistent with the results of Pfenninger and Bacon (1969) but since the 
latter authors gave no details of the size of the disturbances introduced into 
the boundary layer a quantitative comparison between experiment and theory is 
not possible. 

Unfortunately, the expensive nature of the calculation prevented us form 
investigating the effect of suction or blowing on the nonlinearly unstable 
disturbances. It is, of course, possible to investigate particular cases if 
and when experimental results become available. However, it is clear from 
weakly nonlinear theory that the stabilizing influence of suction on 
transition suggested by the results of our linear calculations of I is perhaps 
destroyed by nonlinear effects. We refer to the fact that, although the 
linear critical Reynolds number increases with suction, the part of the 
neutral curve where subcritical disturbances exist increases. Hence, if 
transition is in any way related to the subcritical disturbances, then the 
suction leads to a larger band of nonlinearly unstable modes. In contrast, if 
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th e blowing at the wall is sufficiently large then only supercritically 
bifurcating modes can exist and no nonlinearly unstable modes will exist. Of 
course, this discussion ignores completely the role of three-dimensional 
disturbances in the instability-transition process so that perhaps the most 
that we should assume is that the stabilizing or destabilizing effect of 
suction on the attachment line instability problem is not yet as fully 
understood; obviously the present calculation suggests many experimental 
aspects of the problem which have not yet been investigated. 

Unfortunately, the particular X-dependence of the problem that we assumed 
in Section 2 does not generalize to oblique waves so that this type of 
disturbance can only be investigated in a formally self-consistent manner 
using asymptotic means based on a high Reynolds number assumption. In any 
flow of practical importance the basic flow which we have considered is only 
the first approximation to the flow near the attachment line. It is yet to be 
shown how the Tollmien-Schlichting instability along the attachment line 
merges into a 'crossflow' instability further away from the attachment line. 
There again it seems that the only self-consistent way to investigate this 
problem would be to use asymptotic methods based on a high Reynolds number 


assumption. 
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Table I. Comparison of Calculated Growth Rate with Linear Theory 

(K - 4, t » 306, e = .0001) 


R 

a 

a linear 

Calculated 

| Error | 

570 

.25648 

-.0004523 

-.0004365 

.0000158 

610 

.25648 

0 . 

.0000160 

.0000160 

655 

.25648 

.0004423 

.0004589 

.0000166 


Table II. Summary of Navier— Stokes Computations for Subcritical Instability 

(g = grows, d = decays, n = neutral) 


a 


.28 

.3 

.32 

.325 

.33 

.335 

.34 

.3425 

.345 

.35 

.37 

.4 


R 


570 560 550 540 


d 

g 

g 

g 


g 

g 

g 


g 

g 


d 

n 

S 

g 


g 

d 


535 530 500 


d 


d 

d d 

n d 

d 

d d 

d 
d 
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Figure Captions 

Figure 1. The neutral curves in the ct Q - R Q plane for different values of 

K . 

Figure 2. The neutral curves in the °q c q “ Rq plane for different values 
of K. 

Figure 3. The real parts of the constants a Q and a l for different values 
of k . 

Figure 4. Comparison of the calculated neutral curve with the experiment. The 
finite amplitude solution bifurcates subcritically beyond the arrow 
marked along the upper branch. 

Figure 5. A plot of the (a) eigenfunction, (b) first harmonic and (c) mean 
flow correction corresponding to the linear critical point for 
< = 0 . 

Figure 6. A plot of the calculated wall pressure for R = 570, a = .32, 
e = .00001. Here £ is the initial perturbation amplitude. 

Figure 7. Computed results for R = 570, a = .28, e = .12. 

(a) disturbance energy 

(b) disturbance growth rate 

(c) amplitude of fundamental mode (1) and the first harmonic (2) 

(d) wall pressure 
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Figure 8. Computed results for R = 570, a = .34, e = .12. The legend is the 
same as in Figure 7 . 

Figure 9. A plot of disturbance energy as a function of time for o = .34 
and c = .12 for four different Reynolds numbers. 

Figure 10. Computed results for R = 500, a = .35, e = .2. 

Figure 11. Computed results for R = 900, a = -.201178 and e = .02, .04, .06. 

(a) disturbance energy 

(b) disturbance growth rate 

(c) amplitude of the fundamental mode 

Figure 12. Computed amplitude of the fundamental mode for R = 1000, 
a = .201178 and e = .05, .065. 
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Figure 4 
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Figure 5a 
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Figure 5b 
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Figure 7 a 
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Figure 10a 
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